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Abstract 

The correlations of the fluctuating stress tensor are calculated in an equilibrium molecular-dynamics simulation of a 
Lennard-Jones liquid. We define a coarse-grained local stress tensor which can be calculated numerically and which 
allows for the first time to determine the stress correlation function both in time and in space. Our findings corroborate 
the assumptions made in fluctuating hydrodynamics as long as the liquid is isotropic, that is in bulk. In the vicinity of a 
rigid plate, however, the isotropy is restricted, and major modifications must be done with respect to the usual theory. 
Among these are the appearance of five different viscosities instead of two and a non-trivial dependence of the distance 
from the wall. We determine these viscosities from the simulation data and find that their values are very different from 
the bulk values. We further find much longer relaxation times of the stress correlations than in bulk. 
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1. Introduction 

The theory of fluctuating hydrodynamics starts with 
classical hydrodynamics, which describes the average mo- 
tion of a fluid and which is formulated in terms of smooth 
tensor flelds, and enriches it by additive fluctuations. It 
can be understood as the attempt to bridge the gap be- 
tween microfluidics and nanofluidics from the large-scales 
side. The central question is of course, what properties 
the fluctuations must have in order to correctly represent 
the microscopic states of the liquid in a given problem? 
' The standard approach [1] [2] is to split the full veloc- 
ity and stress tensor into a macroscopic (denoted by an 
overbar in the following) and into fluctuating parts (de- 
noted by a tilde). The linearised governing equations take 
then the form of a two-scale process, with the standard 
compressible Navier-Stokes equations for the macroscopic 
velocity v"(x,i) and the corresponding Newtonian stress 



tensor cTij (x, t) , 



Here and in the following, 77 denotes the shear viscosity of 
the fluid and A its volume viscosity. As we will be con- 
cerned only about equilibrium situations, the macroscopic 
equations reduce to 



p(x,t) 



p(x,i) 



(2) 



po, pix,tj=po, v(x,t) = 0, 

with constants po ^nd po- The fluctuating flelds satisfy a 
linearised Navier-Stokes equation with additional random 
noise sa, 



Podtv = div(a--l-s), 
dtjj = -podivv. 



(3) 
(4) 



The stress tensor in equation ([3]) is again split into two 
parts, the flrst of which satisfies the same Newtonian def- 
inition (fTl) as the macroscopic stress tensor, only in terms 
of a fluctuating velocity v(x, t) and of a fluctuating pres- 
sure p(x,i). The second contribution is an uncorrelated 
fluctuating stress tensor Sij (x, t) . The tensor aij is meant 
as a coarse-grained stress which is fluctuating, but not as 
randomly as Sij, see for example the book by Kubo et al. 
[3] for coarse-graining concepts, and Ref. [1] for substi- 
tute Markov processes. The same two-step separation has 
been used by Hauge and Martin-L6f [5^, Chow and Her- 
mans [B], but without commenting on why it is necessary. 
Equation ([3| has the structure of a Langevin equation, 
with Sij the random driving. This random stress tensor 
vanishes on average, (s,y(x, t)) — 0, and it has to satisfy 
the fluctuations-dissipation theorem of the second kind. It 
is thus uncorrelated in space and in time, 

(s,j(x,i)sfc,(y,s)) = 2kTA,^kiSi^-y)6{t-s). (5) 



'Stjp{x,t) + r]{diVj+djVi) + (\--rASij divv. (1) Here, (•) 
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Here, (•) denotes the average over different realisations, or 
a time average; k is Boltzmann's constant, and T the tem- 
perature of the system. The correlator is further assumed 
to be isotropic, which flxes the fourth-rank tensor Aijki 

to [a 

A-ijki = vi^ikSji + SiiSjk] + [^- i^mSijSki- (6) 

This special form of Aijki follows from the fact that it 
must be isotropic, thus invariant under any rotation of the 
coordinate system. It is therefore a linear combination of 
the only three isotropic tensors of rank four, which are 
SijSki, SikSji, and SuSjk [Zj- The factors of proportion- 
ality are found to be the three viscosities of a fluid [8J. 
The additional symmetry of the stress tensor, leading to 
Aijki = Ajiki — Aijik eliminates the rotational viscosity. 
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In the vicinity of a rigid boundary, however, isotropy is 
not granted in the same way as in the bulk. Instead, the 
stress correlator should reflect the fact that the normal 
direction is different from the two tangential directions. 
In this paper, we therefore ask which tcnsorial form the 
correlator will adopt when evaluated close to a rigid wall. 
In order to start with the most simple situation, we regard 
the stress correlator in thermodynamic equilibrium. 

We must further ask, what "close" means, that is up 
to which distance from the wall the different symmetry is 
sensed by the fluctuating part of the stress tensor. Near 
the boundary, we thus question also the spatial delta func- 
tion in Eq. (Is]). The problem with the spatial delta func- 
tion becomes apparent when applying fluctuating hydro- 
dynamics for example to the Brownian motion of a spher- 
ical particle. This has been done in the literature on the 
autocorrelation of a Brownian particle [51 15] . If we naively 
integrate the random stress tensor over the particle in or- 
der to obtain the random force, and then construct the 
autocorrelation of that force, the result would always be 
infinite. The reason is that we then integrate the three- 
dimensional delta function from Eq. ([5]) only over the two- 
dimensional surface of the particle. The third dimension 
thus remains a delta function, evaluated at zero. 

The spatial delta function must be questioned also from 
a less formal point of view: The standard argument for the 
absence of spatial correlations is based on thermodynam- 
ics [T] : The fluid volume is cut into several small but finite 
portions AV. In each of these volumes, the entropy pro- 
duction is calculated as the volume-integral of the tensorial 
reduction of the stress with the rate-of-strain. The fluctua- 
tions of the stress tensor is then taken from the variations 
of the entropy production around the equilibrium state. 
As this integral does not contain any surface contribution, 
the fluctuations of two stress tensors attributed to two such 
boxes are independent of each other. This property is then 
kept even in the limit of infinitesimal box size. It appears 
as if there should be a lower spatial bound to the validity 
of this argument, namely when the volume size become of 
the order of the distance between the molecules of the liq- 
uid. As thermodynamic arguments are used, the thermo- 
dynamic limit of many interacting molecules is assumed to 
be done before the box size can be sent to zero. It appears 
thus as a reasonable question to ask, down to which spatial 
scale does the delta function yield a good description? It 
will turn out in the following, that it works amazingly well: 
In bulk fluid, the delta function holds down to distances 
smaller than the intermolecular spacing. 

It is the aim of this paper to test the validity of Eqs. ^ 
and ([6]), describing stress correlations — in particular, to 
find deviations in the vicinity of a rigid surface. In or- 
der to do so, we start from a fully microscopic simula- 
tion of a liquid and extract the stress tensor using coarse- 
graining. As the simulation is done at thermal equilib- 
rium, the macroscopic fluid velocity v vanishes, leaving 
only the thermodynamic pressure in the macroscopic stress 
tensor aij — —poSij. We correlate the coarse-grained stress 



tensor at different positions and at different times in or- 
der to verify both the isotropy in Eq. ^ and the delta- 
functions in Eq. ([5]). This will be done both near a rigid 
wall and in bulk ffuid. 

The precise form of the stress tensor correlations is po- 
tentially important for all ffuidic systems including bound- 
aries. In microfluidic applications, and even more in nano- 
fluidic ones, the role of the boundaries is essential [TOl 
im 112] . In nanofluidics, the presence of boundaries even 
introduces different time scales in the relaxation dynam- 
ics [13]. For example, Brownian motion of a particle im- 
mersed in a viscous fluid, which is one of the most classi- 
cal systems both in the hydrodynamic and in the stochas- 
tic community, is in most cases a microfluidic problem — 
depending on the size of the particle. Both for large and 
for small particles, the velocity fluctuations are known to 
exhibit long-time tails in their autocorrelation-function. 
The theoretical explanations of this phenomenon, as far 
as they consider both the particle and the fluid, assume 
the isotropic form of Eqs. ^ and (|6| for the fluctuating 
part of the stress tensor — even when this stress tensor is 
integrated over the surface of the particle to yield the total 
fluctuating force acting on it [3] |B] . The investigation in 
the present paper is one of two necessary steps to render 
the theoretical description of the long-time tails more con- 
sistent. We will flnd below that it is not the usual viscosity 
which occurs in the formula of the force autocorrelation, 
but a viscosity which is a property of the fluid and the 
structuring effect of the rigid particle itself. The other 
step will be published in Ref. [14j . 

We will further flnd that the inclusion of boundaries 
makes it necessary to pass from homogeneous viscosities 
(which are deflned in the limit of homogeneous space with- 
out any scale) to rather local viscosities. Such viscosities 
reflect more details of the geometrical conditions and of 
the time scales involved. 



2. Symmetry considerations 

The form ^ of the stress correlator follows directly 
from symmetry considerations. We will now apply the 
same considerations to derive the corresponding expression 
for a liquid close to a wall. Similar formulae are found 
in the literature on nematic liquids |15| . The following 
symmetries apply to the stress correlator: 

(A) The most rigorously enforced symmetry is that the 
stress tensor is sym,metric. This symmetry is inherited 
from the microscopic deflnition of the stress tensor and is 
thus not subject to random fluctuations. It leads to the 
conditions 



(SijSkl) = (SjiSkl) = {SijSlk). 



(7) 



By this symmetry, the 81 components in the full corre- 
lation tensor of rank four in three-dimensional space are 
reduced to 36 different ones. 
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Figure 1: Schematic breaking of viscosities in bulk fluid (left) into 
the surface viscosities near a flat wall (right). Given are the three/six 
non-zero values found in the components of the stress tensor correla- 
tor Aijki, indexed by the normal direction N and the two tangential 
directions T and T' . 



(B) Another symmetry is isotropy, that is the invari- 
ance of the correlator under rotations. It is not enforced 
by the microscopic definition and must result from the av- 
eraging either over an ensemble or possibly over a long 
trajectory. In bulk, the rotations are arbitrary, while in 
a geometry with plates one direction is distinct from the 
two others and only rotations parallel to the plates are per- 
mitted. The maximal number of different components is 
now reduced to 21. They are, however, not all completely 
independent, and some of them vanish: A systematic de- 
composition of the rank-four tensor into parallel and or- 
thogonal parts, using parallel isotropy and Eqs. ([t]), reveals 
that these 21 entries can take only nine different non-zero 
values, and that they depend on seven independent pa- 
rameters (viscosities) . 

(C) The next symmetry is the invariance under index 
pair exchange^ 



(SijSkl) = (sklSij). 



(8) 



It is a consequence of the entropy production being a 
scalar [8]. This symmetry leads to the identification of 
two of the viscosities and the vanishing of another. After 
having used all symmetry requirements, we thus have a 
total of five viscosities left, which give rise to six different 
non-zero values in the 13 principally different components 
of the stress correlator. 

The symmetry conditions near a plate require remark- 
ably more viscosities than those in bulk, where two vis- 
cosities lead to three difi^erent values. The presence of the 
wall thus induces a separation of the bulk viscosities into 
surface viscosities. This splitting is depicted in Fig. [I] The 
diagram provides on the left-hand side the three (in bulk) 
different values found in the components of Eq. ([6]), and 
on the right-hand side the six non-zero values near a fiat 
rigid immobile plane. The full rank- four tensor expressing 
the random stress correlations near a wall is 

Ajki = Vi NiNjNkNi + ?72 gij9ki 
+ m {9ikgji+gii9jk-gij9ki) + m (NiNjgki+NkNig^j) 
+ V5 iN^Nkgji+NjNkga+N,Nigk+NjNig,k), (9) 



The rightmost column of Fig. [T] uses a symbolic nota- 
tion for the indices of tensor components. A^ stands for 
the direction normal to the wall, T may stand for both 
tangent directions, and T" is only used if both tangent 
directions are used: In particular, TT stands for the two 
index combinations with the same tangent vector, whereas 
TT' stands for the ones with different tangent vectors. 
The same notation will be used in the figures of Sec. |4] 
to visualize the 13 components of the stress correlation 
tensor, six of which are non-zero and seven vanish. 

3. Coarse-graining of a microscopic stress tensor 

Let us consider a molecular dynamics simulation in or- 
der to study the microscopic dynamics of the molecules in 
a liquid. The liquid here consists of point-like atoms, the 
dynamics of which is governed by Hamilton's equations 
comprising a pair-interaction potential. In this dynamics 
we define a microscopic stress tensor. It is chosen such 
that its coarse-grained version, which is obtained by in- 
tegrating it over a small region, will be consistent with a 
macroscopic notion of the hydrodynamic stress tensor. 

The Hamiltonian of N point-like atoms (a) of masses 
m^°'\ positions r(")(i) and linear momenta p(")(i), which 
undergo mutual infiuence via a pair-potential </>, reads 

^ (a) (a) ^ "-1 

^ = E^^;;^ + EE^(l|r^'^-^"'ll)- (10) 



a=l p=l 



The time evolutions of the positions and momenta are 
given by Hamilton's equations. The microscopic stress ten- 
sor Tik, can then be introduced as the current density of 
the density of linear momentum g, reading 



5.(r,i):=^^pf^(t)<5(r-r(")(t)). 



(11) 



In the absence of external forces, the conservation of mo- 
mentum is expressed by the following balance equation, 
which introduces the microscopic stress tensor r^fc. 



dtgi{-r,t) -dr^nk{r,t) = 0. 



(12) 



The definition of these tensor fields allows us to pass from 
a description in terms of the atom positions r(")(i) to a 
field description for any position r. Equation ( [T2| , which is 
at the same time a conservation theorem and a definition, 
fixes r only partially. We may always add an arbitrary part 
with zero divergence, a gauge freedom which is inherent to 
the continuum description. However, we will not go into 
the implications of this freedom. A convenient definition 
of the microscopic stress tensor is the following, since it 
is symmetric and lends itself to a straight-forward coarse- 



with the shortcut gij 



N,N, 



graining technique, 



m 



'EE^'di^^'^-^" 



a P^a 



l^||r(/3)-r(")| 



(^M.^t/S)) /"5(r-r(")-A(r(«-r(")))dA. (13) 



The first term is usuaUy caUed the kinetic part of the stress 
tensor. The second term comprises virial contributions, to- 
gether with the integral of a three-dimensional delta func- 
tion over the line starting at the point r^") and ending 
at r*-^-'. Expression (13 1 is the same as that of Forster [HI 



see Eq. (4.6)], with the exception of the lower boundary of 
this integral. 

We choose to do a coarse-graining by averaging the 
microscopic stress tensor within a small cube y(x) of 
side-length a, centred around x. We call this average the 
coarse-grained stress tensor t^^ (x, t) at the point x. 



'ik 



(x,i) 



1 



Tik{r,t)dr. 



(14) 



y-Cx) 



The integration of the kinetic term in Tik is readily done. It 
uses the momenta of the particles which are found within 
the volume at time t. The line integral of the delta function 
has a simple geometrical interpretation. When integrated 
over the volume V"^ (x) , only those parts of the line integral 
contribute where the line is found within the volume. For 
the (convex) cubic box y (x), only four different cases are 
possible: 



||i.(/3) _!.(") II / / ,5('r-r(")-A(r('3)-r("))) dAdr 



V 

|r(/9)_r(")| 
l|r-r(")|| 



if r(°) G V and r^'') e V, 
if r(") e V^ and r^'') ^ F, 
||r(«_f II if r(") ^ V and r(^) e V, 
||r-f'|| if r(") ^ V^ and r(^) ^ F. 



(15) 



Here, f and r' denote the points where the connection 
line between the atoms intersect the boundary of the vol- 
ume y°(x). The last case in Eq. (15) yields no contri- 



bution if these intersection points do not exist. Figure [2] 
provides a graphical representation of these different cases. 
For calculating the coarse-grained stress tensor we thus 
have to find the contributing atom pairs, to evaluate the 
pair potential, and to find the intersection points. Notice 
that even a box which does not contain any particle may 
yield a non-zero microscopic stress. It is sufficient that any 
connecting line passes through the box. 

The usual way to calculate a stress tensor in molecu- 
lar dynamics simulations is to add up the kinetic contri- 
bution and the full virial of the atoms. The simulation 




Figure 2: Visualisation of the virial contributions to the coarse- 
graining of the stress tensor. Prom the connection lines between 
the atoms (small circles) only the bold parts are taken into account. 
All cases of Eq. ||15|l are displayed. 



package LAMMPS [T7|, for example, outputs the following 
quantity as stress per atomPJ 



T-(a) (a) (a) (a) , 1 V^ / (a) (0)\ 

Vi. ■■= -m'^ 'v) vl ' + - 2_^ (r^ -^r ) ^ 

;geNbrs.(Q) 



■"l "k 



i i 

||r(«)_r(/3)|| 



,'nir(°)_r(^)| 



0'(l 



(16) 



-(a) 



However, T^f. has the dimension of stress times volume. 
The stress is attributed to the atom with index a, its neigh- 
bours have indices (3. It is clear that 7^^" cannot be a 
local quantity, as it does not depend on any space vari- 
able. It may be the integral of a local quantity — such as 
the stress in Eq. ( 13 ) — over a specific volume around the 
atom a — but over which volume? We find the answer to 
this question from the explicit averaging procedure above: 
Only if the particle and all its neighbours are within the 



volume of integration, that is the first case in Eq. (15), 



then we obtain expression (16) as the contribution com 



-(a) 



ing from this atom. We conclude that 7J^ ' is the bulk 
contribution to a stress tensor when integrated over a vol- 
ume which is sufficiently large that it contains the atom 
in question and all its interaction partners. It discards all 
boundary terms. For geometrical reasons it is not possi- 
ble to assign space-filling individual volumes to all atoms 
(such as is done for example by a Voronoi tessellation) in 
order to obtain the stress per atom as an integral over the 
local stress. As the cell boundaries in the Voronoi tessella- 
tion are determined only by nearest neighbour's positions, 
they cannot correctly represent all interaction partner's 
positions, which are spanned over a larger region. 

As our interest lies in the local fluctuations of the stress 
tensor, we have to keep the integration boxes as small as 



possible. The small boxes forbid the use of Eq. (16), we 



have to average the microscopic stress tensor explicitly. 



according to the cases in Eq. ( 15 ). The additional work to 



be done in comparison with Eq. ( 16 ) is the solution of the 
geometrical intersection problem. 

^Throughout the paper, terms specific to LAMMPS are typeset 
in a typewriter font. 




(b) 
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Figure 3: (Colour on-line) Snapshots of the simulation domain with 
N = 1380 fluid atoms (black), confined between two fixed plates 
which consist of a regular grid of immobile atoms (red/blue). Shown 
are cuts through the domain as seen from two different sides. The 
series of cubes in the middle of panel (a) indicate the cubic boxes 
used for averaging in Sec. |4.2| In the hatched region, a Langevin 
thermostat controls the temperature. 



4. Results from molecular dynamics simulations 

The molecular dynamics simulations described in the 
following have been done with the molecular dynamics 
package LAMMPS il7:. The coarse-graining procedure de- 
scribed above has been added as a module. 

4.I. Details of the simulation 

Our simulation is set up as a three-dimensional Lennard- 
Jones liquid between two plates, as depicted in Fig. [3] 
1380 fluid atoms are confined between two fixed plates. 
The plates are constructed as a single regular layer of 
atoms with the same interaction potential. In order to 
approximate a smooth flat surface, the lateral distance 
between these plate atoms is about a third of that in the 
fluid. The simulation domain is periodically continued in 
the two directions parallel to the plates. 

The fluid atoms are advanced in time using the Verlet 
algorithm (NVE integration) which is consistent with the 
microcanonical ensemble. An exception are the hatched 
regions in Fig. [Sk, where an additional Langevin thermo- 
stat is applied. Such a thermostat is necessary to counter- 
act numerical integration errors, but it modifies the statis- 
tics of the stress. We therefore chose to apply it only in 
regions where we do not measure the stress. The plate 



atoms are not advanced in time, simply by not applying 
any integration routine. The forces on these atoms are 
calculated, however, which is equivalent to the possible in- 
clusion of strong constraining forces from crystalline layers 
behind. This procedure allows to determine the force ex- 
erted on the plates by the fluid. 

The pair-interaction potential between the atoms is the 
12/6 Lennard- Jones potential, modified by an interpola- 
tion term to have a finite carrier. 



0(r) 



4e[{a/ry'-ia/rf] 




for r < 2.5a 
for r > 3. 5(7 

else. 



(17) 



n— — 3 



The six parameters ^„ are determined to guarantee con- 
tinuous derivatives everywhere, up to the second. We use 
this smooth potential in order to make sure that discon- 
tinuities in the second derivative do not produce artefacts 
in the measurement of the stress tensor. 

In the following, we will use Lennard-Jones units for 
all quantities. In these units, distances are measured as 
multiples of the length parameter a in the Lennard-Jones 
potential, energy as multiples of e, and time in terms of 
•^TOcr^/e, where m is the mass unit. In these units, the 
mass of the atoms is chosen to be 1, the side-length of 
the simulation domain is fixed to be 12. The time unit 
is approximately the oscillation period in the minimum of 
the pair potential (/). The time-step of integration, used in 
the Verlet algorithm is 0.002, thus more than two orders 
of magnitude smaller than all physical time scales in the 
system. 

The density and the (kinetic) temperature of the sys- 
tem are chosen such that it is found in its liquid state, 
according to the phase diagram in Refs. [HI [191 120]. Of 
course, this holds only in the "bulk" region where the 
structuring effect of the walls is absent. The disordered 
fluid is visible in the middle of Fig. [3^. There, the follow- 
ing average values for temperature, density and pressure 
are found: T — 1.04, po = 0.83, po = 1-3- Near the plates, 
the enforced flat geometry imposes a partial ordering of 
the fluid. Up to three distinct layers of fluid atoms can 
be identifled near each wall. The first layer has a locally 
hexagonal structure as seen in Fig. [Sb. Its width and ori- 
entation are independent of the grid spacing of the plate 
atoms. Due to periodic boundary conditions, the hexag- 
onal grid does not rotate anymore as soon as the system 
is equilibrated. Of course, also the stress tensor becomes 
anisotropic close to the plates. 

The simulation was first run to let the system find its 
thermodynamic equilibrium. We then continued the simu- 
lation while calculating the coarse-grained stress tensors at 
some fixed positions every 5 time steps. These time series 
are analysed in the following, either using the blocking 
method |5T] to ensure that all used variables are equili- 
brated, or by simple averaging to measure spatial correla- 
tions of the stress tensor, or by fast Fourier transforms for 
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thus find the deviations of the stress tensor to be 



a in LJ-units 

Figure 4: (Colour on-line) Fluctuations of the stress tensor in bulk, as 
a function of the side a of the averaging box. All tensor components 
are plotted. They are grouped according to the expected symmetry 
(see text). 




number of data points in the "blocked" list 

Figure 5: (Colour on-line) Estimated error of the average estimation 
from the blocking transformation |21| . We used the time series of 
all stress components in a box of size a = 1.0 centred between the 
plates. 



the time correlations. The averages (•) in the following are 
thus time-averages over a long trajectory. As the system 
is in equilibrium, this should be equivalent to an ensemble 
average. 

^.2. Varying the size of the coarse-graining boxes 

Before presenting the temporal and spatial dependence 
of the stress correlator, let us investigate how the side 
length a of the averaging box changes the fluctuations of 
the coarse-grained stress tensor inside. 

The explicit coarse-graining method by averaging over 
small cubes allows to observe the transition from thermo- 
dynamic behaviour (in large boxes) to highly fluctuating 
random behaviour (in small boxes). The amount of fluc- 
tuations depends of course on the size of the box. As we 
treat here only thermodynamically equilibrated fluids, the 
mean square deviation of the stress tensor should scale as 
the inverse volume of the averaging box - just as for any 
intensive quantity ;22, § 112]. In the bulk fluid, we should 



\/(^JW)^««"'^'' 



(18) 



where fj " (x, t) is defined by subtracting the macroscopic 
stress tensor, which is —poSij, 



rf)(x,i):=rJ^(x,t)+po'5, 



(19) 



The value for the pressure is taken from disordered fluid 
atoms only, and not from the layers close to the plates. 



The numerical verification of the scaling ( 18 ) is plotted 



in Fig.|4] Amazingly, we find the correct scaling already for 
box sizes of a ^ 2, that is a side length of the order of the 



average distance between the fluid atoms. The scaling ( 18 1 



is originally based on a thermodynamic argument which 
requires that many particles are involved. Instead, we here 
see the correct scaling already for very few particles. It 
appears that there is more truth in the thermodynamic 
argument than just the interaction of many particles. 

For visualising the tensor of stress fluctuations, we have 
chosen the following strategy, which will be applied also to 
the following plots: The symmetry requirements (B) and 
(C) in Sec.l2]lead to components which are equal only in a 
statistical sense. All components which may exhibit such 
fluctuations are plotted by individual lines. The legend 
and the colour coding, however, differentiate only between 
the groups of components which should be equal on av- 
erage. The number of lines actually plotted is given in 
parentheses. 

The fact that there is no difference between the nor- 
mal and the parallel directions in the data of Fig. |4] con- 
firms our assumption that the fluid behaviour in the mid- 
dle between the walls is the same as if they were absent. 
There, the diagonal and the off-diagonal parts fluctuate 
differently. This is a consequence of the components being 
different functions of the two viscosities (see also Fig. IT]). 

Exemplary for all plotted variables, we show for the 
centred boxes that the system is well equilibrated. Fig. [5] 
visualises the blocking transformation [3T] of the diagonal 
entries in the coarse-grained stress tensor. A plateau is 
clearly visible, which shows that the analysed time series 
indeed correspond to a stationary stochastic process. 

4-3. Stress correlations in bulk fluid 

The aim of this paper is to determine the spatial and 
the temporal structure in correlations of fluctuations of 
the coarse-grained stress tensor f^^ . We are interested in 
the difference between correlations in bulk and in the near 
vicinity of the walls. In order to do so, we set out three 
different series of cubic boxes, as depicted in Fig.[6j These 
cubes are shifted by small distances and may overlap. The 
first series is in the middle between the plates, where the 
fiuid behaves as in bulk, a second one close to a wall and 
parallel to it, and a third series is normal to the wall. 

The displacement of the boxes makes it possible to 
determine the spatial form of the stress correlation by 
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Figure 6: (Colour on-line) Same as Fig. pa, but with different series 
of overlapping coarse-graining cubes, which are used in Figs. [TjfTo] 



cross-correlating the coarse-grained stress fluctuations in 
different boxes. Temporal correlation functions will also 
be given. We would like to compare the result with the 
uncorrelated theory in Eq. ([5|. However, the nonzero size 
of the boxes used for averaging does not allow a direct 
comparison. Remind that fj" depends on the box size. 
We should find the delta function of Eq. (Is]) in the limit 
a — > 0, which is infeasible in practice. For all nonzero 
sizes, the numerical data is sort of "convoluted" with the 
shape of the boxes. The best we can do is to compare with 
the equally smoothed version of Eq. ^ using a reasonably 
small a. The corresponding correlator of the uncorrelated 
stress tensor Sy is expressed by the relative overlap of the 
two boxes. 



{r'^frif)^\V%^)lJV%y% 



(20) 



of course disregarding the infinite term d{t—s) in Eq. ([5]). 
Here, | • | denotes taking the volume of a set. 

The spatial and the temporal correlations of the ran- 
dom coarse-grained stress tensor are depicted in Fig. [7] 
The data have been generated by correlating the contents 
each box in the central series of Fig. [6] with the first one, 
which has Xq as its centre. Plotted are all 36 components 
of the resulting rank-four tensor which can be different us- 
ing the symmetry (A) of Sec. ^ They are grouped into 
the 13 different classes according to the symmetries (B) 
and (C). Fig. [T] shows that seven of them vanish (grey 
lines) and that the remaining six take three different val- 
ues. This is expected according to what has been said 
around Fig. [l] 

The spatial correlations in Fig. [7^ should be compared 
with the results expected from an uncorrelated stress ten- 
sor. These would be lines starting at some value which 
is related to the viscosities (unspecified due to 5{t—s) in 
Eq. ([5])) and which goes to zero the distance a and beyond. 
The coincidence in Fig. [Tk is remarkable. The numeric 
values are a bit more smoothed, but the apparent convo- 
lution kernel responsible for this smoothing has a width of 
less than unity. We therefore conclude that Eqs. ([5]) with 
Eq. m in bulk are very precise, down to a length scale of 



the order of the molecule distance. Similarly to what we 
found in Sec. |4.2[ such a strong result cannot be expected 
from thermodynamic considerations. 

The temporal correlations in Fig. [7)d exhibit the devi- 
ations from the temporal delta function in Eq. (Is]). The 
decay is sufficiently fast to reasonably adopt the approx- 
imation of the delta function. Whether there is an alge- 
braic contribution, the so-called long-time tail cannot be 
judged from these data. In order to calculate the viscosity 
values, we use the integral over the detailed time correla- 
tor, as proposed by the Green-Kubo equation. The idea 
is the same as in Refs. [2311211 Hi] j only that we take a box 
smaller than the simulation box. For the shear viscosity, 
this reads 



kT 



dt(fi,"^(x,o)f(;^(x,f) 



(21) 



and similarly for the two components comprising the vol- 
ume viscosity A. From the data in Fig. [Tb we find the 
values rj = 0.84 and A = 1.72 by regression^ 

4.4- Stress correlations near a rigid flat wall 

Near the plates, the stress correlations look quantita- 
tively and qualitatively different. The average stress ten- 
sor and also its fluctuations are anisotropic. Our aim in 
this section is to determine the five surface viscosities in 
order to quantify their deviations from the bulk values, 
according to the scheme in Fig. [T] 

First, let us look at the coarse-grained random stress 
tensor itself. Fig. [8] reveals the mechanical properties of 
the first fluid layers: We find a huge lateral diagonal stress 
( TT components) which strongly depends on the distance 
to the wall but not on the lateral position. In the fig- 
ure, the points x correspond to the centres of the series 
of boxes laid out in Fig. [6j The position Xq is the centre 
of the first box in the normal series, which is closest to 
the wall. The average value of the huge diagonal stress is 
even larger than its fluctuations. This lateral stress can 
be understood as a surface tension of the contact between 
liquid and solid. The AW-component is much smaller — its 
average value is the same as the bulk pressure. The A^A^- 
component is dominated by its fluctuations which grow in 
the vicinity of the wall and even overgrow the tangential 
stresses. In the following, we will flnd that there are cor- 
relations in these fluctuations. The lines in Fig. [Sjs prove 
that both the average values and the standard deviations 
are independent of the lateral position. The roughness of 
the AW-component in this plot is due to the large fluctu- 
ations. 



■^Notice that these values do not quite correspond to the finding 
in Refs. |24l 125) . A possible explanation is the dependence of the 
viscosities on the size a of the averaging box. We found such a 
dependence together with long-time tails in the autocorrelation of 
the random stress tensor in bulk. 



7 




NN-NN{lx) 
TT-TT{2x) 
TT-^T'T' (2x) 
NN-TT(4x) 
TT'-TT' (Ix) 
NT-NT i2x) 
NN-NT{4x) 
TT-TN (Ax) 
TT-TT' {4x) 
NN-TT' (Ix) 
TT-T'N i4x) 
NT-NT' (2x) 
TN-TT' {4x) 



NN-NN{\x) 
TT-TT(2x) 
TT-T'T' (2x) 
NN-TT (4x) 
TT'-TT' (Ix) 
NT-NT (2x) 
NN-NT(4x) 
TT-TN (4x) 
TT-TT' (4x) 
NN-TT' (2x) 
TT-T'N (4x) 
NT-NT' {2x) 
TN-TT' (4x) 



0.4 0.5 0.6 

|t — s| itiLJ-units 



Figure 7: (Colour on-line) Spatial and temporal correlations of the stress fluctuations in the middle of the domain (in bulk). Plotted are the 
average correlators and their standard deviations (S(-), see insets) for different classes of index combinations. 



For the data in Fig.|8]and in the foUowing ones, another 
comment on the symmetry is necessary. By the noted sym- 
metries in Sec. [21 the two hnes for the TT-components in 
Fig. [8] should coincide. However, they do not. The reason 
is that the structured hquid is not entirely isotropic in the 
tangential directions, but it is ordered in a hexagonal grid, 
as depicted in Fig. [Sjs. This strict ordering is a result of 
the wall being very flat. A roughness of the size of one 
atom distance would easily destroy the hexagonal struc- 
ture, and only the symmetries of Sec. [2] would remain. In 
order to weaken the effect, we averaged over eight trajecto- 
ries which differed only by the orientation of the hexagonal 
grid. This reduced the deviations considerably to those in 
Fig.i 

The spatial correlations of the random stress tensor 
are depicted in Fig. [9) Again, the two series of integra- 
tion boxes in different directions are used, see Fig. |6] The 
plot shows the correlation of each cube in the series with 
the first one, which is in the normal series the one closely 
aligned at the wall. For the normal series, the result is not 
expected to be independent of this choice. Only in the par- 
allel series, the correlation should depend only on the dis- 
tance of the boxes. As explained above, there are in prin- 
ciple 13 different results to be expected. Seven of them, 
those indicated in grey, are zero. The six non-zero compo- 
nents are plotted in colour. The correlations of the stress 
tensor inherit the main features of the stress in Fig. [81 in 
particular the systematic part which we called a surface 



tension. In normal direction, the tangential components 
exhibit a similar anti-correlation as in Fig. [8}d. The decay 
of the most pronounced correlations takes place on the 
length of one atom-distance and is independent of the box 
size a. We verified this statement using other box sizes. In 
tangential direction, the surface tension manifests itself in 
a nonzero value of the tangential correlation, see Fig. [8b. 
Apart from the surface tension, we find in Fig. [8p a 
similar behaviour as in bulk, that is a linear decay up to 
a distance a and zero beyond. This is clearly visible for 
the NN-NN component, but also for the others — after 
subtracting the contribution of the surface tension. This 
implies that the spatial delta function still works parallel 
to the walls. The correlation in normal direction, which is 
dominated by the details of the interaction between wall 
and fluid, depends both on the distance between the box 
and their distance to the wall. These two considerations 
lead us to the following proposition for the correlation 
function of the random stress tensor near an immobile flat 
rigid wall: 

(s.y(x,t)sfe;(y, s)) = 2kT Aijki5{t-s) x 

<52(xll-yll)/(x^,y^). (22) 

with the tensor Aij^i from Eq. ([9|), depending on the sur- 
face viscosities, and with x" and x-'- the parallel and nor- 
mal projections, respectively. The function /, which has 
the physical dimension of inverse length, takes into ac- 
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Figure 8: (Colour on-line) The coarse-grained stress tensor as a func- 
tion of the position, as indicated in two series close to the wall in 
Fig. p] (a) parallel to the plate, (b) orthogonal to the plate. 



count the microscopic details of the interaction between 
wah and fluid. As the only length scale involved is the 
distance between wall and the first structured liquid layer, 
we expected / to scale as one over this distance. 

It remains to check whether the delta function in time 
is still a good approximation, and what values the surface 
viscosities have. The temporal correlations near the wall 
are depicted in Fig. [TO] The results are taken from one of 
the boxes in the parallel series. Here, we see even better 
than in Fig. [9] that there is a systematic contribution com- 
ing from the surface tension. Notice that it affects only 
the tangential components and not the heavily oscillating 
NN~NN component, which is decaying to zero, but on a 
longer time scale than plotted. 

In order to obtain the viscosities from these data, we 
have subtracted the asymptotic value from the correla- 
tion function, which renders the function integrable. This 
means that we implicitly subtracted a sort of local pressure 
which is allowed to be anisotropic. This pressure reflects 
the static contributions from the interaction liquid-wall, 
which we called surface tension above. The Green-Kubo 
integral of the such normalized correlation functions yields 



the following values for the surface viscosities: 

m « 113, 

r/2 between 16 and 34, 

773 « 8.4, 

774 between 65 and 178, 

775 between 4.2 and 5.5. 



(23) 



The values for these viscosities are not very precise for 
three reasons: The first is the hexagonal structure of the 
first fluid layer which does not coincide with the Cartesian 
coordinates used for the tensor indices. We thus expect 
different viscosity values for different orientations of the 
hexagonal grid. The above mentioned average over eight 
different trajectories with different orientations apparently 
does not suffice to obtain unified values. The second rea- 
son is that the thermodynamic pressure is not isotropic 
anymore and that the time-correlation function does not 
decay to zero. We thus had to strip off a linear time- 
dependency before integrating, as already explained. The 
third reason is that the intrinsic time scales of the time- 
correlation function are much longer than in bulk. While 
the bulk viscosities could be easily identified from a Green- 
Kubo integral over a few dozens of timesteps, it required 
thousand and more timesteps to do the same for the sur- 
face viscosities. This much longer time scale can already 
be guessed from the comparison of Figs. [7)d andfTO] 

4-5. Hydrophobic walls 

The rigid plates used in the above simulations consist of 
the same type of atoms as the liquid. Such walls are known 
to be hydrophilic, which is also seen in the strict ordering of 
the first structured fluid layer. We repeated the simulation 
using hydrophobic walls, which have a modified Lennard- 
Jones potential to interact with the fluid atoms 126) . 



(r)=4e[(a/r)i2-c(a/r)^ 



(24) 



with the same cutoff and with an equivalent smoothing 



of the potential as in Eq. (17). According to Barrat and 



Bocquet i26j, who proposed the form (24) for simulating 
hydrophobic walls, the factor factor c = 0.5 which we used, 
should correspond to a contact angle of about 140°. 

The simulation leads to a similarly structured liquid as 
in Fig. ([3]), only that the first layers are a bit farther away 
from the walls, that they are not as flat, and that their lat- 
eral structure is not as regular. The ordering effect of the 
first layers are nevertheless very strong, and must there- 
fore be induced by the flat geometry of the wall. We took 
the same bulk pressure and density as in the previous sim- 
ulation. As expected for a hydrophobic surface, the TT- 
components in Fig. |8] are different. They are considerably 
smaller, around the value 3 instead of 13. The correlation 
functions look similar to those depicted in Figs. [9] and 10 
with the difference that they decay much faster in time, 
in particular the NN-NN component. Accordingly, also 
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Figure 9: (Colour on-line) Spatial correlations of the random stress tensor in the vicinity of a wall, (a) in normal direction, and (b) parallel 
to the wall, see the series of cubes in Fig. [m Each box in the series is correlated with the first one, centred around xq (in the orthogonal case; 
the one closest to the wall) . Legend and inset as in Fig. [Te.. 



the surface viscosities turn out to be different. Near the 
hydrophobic surfaces, we find the values 

m ~ 87, 
m « 1.7, 

m ~ 1-5, (25) 

r/4 « 8.1, 
rys w 0.35. 

It is not astonishing that these viscosities are smaller than 
those found near the wetting walls. In particular in tan- 
gential direction, where the liquid is feels less constraints, 
the dissipation is smaller. This affects 772~?75, but not 
771, which is accordingly of a similar magnitude as in the 
wetting- wall case. Particularly astonishing is the small 
value of 775, which takes over the role of the shear viscos- 
ity. It is much smaller than the bulk value (r/ = 0.84 also 
in the hydrophobic simulation). 

5. Discussion 

The numerical observations above, as described in Sec- 
tions |43] and |44] allow to estimate the validity of the main 
assumptions leading to fluctuating hydrodynamics. In- 
deed, we observed that the spatial delta function in Eq. ^ 



works very well in bulk. Near a plate, which has been 
assumed to be flat, immobile and rigid, the isotropy is 
restricted to rotations parallel to the wall, and we must 
replace the spatial delta function by one which is only act- 
ing parallel to the wall, as given in Eq. ( 22 ) . In orthogonal 



direction, the details of the dynamical interaction between 
liquid and wall introduce an unknown factor in form of 
a function depending on the orthogonal distances. More- 
over, we have to deal with five viscosities instead of two, 
the values of which are provided in Eq. ( 23 1 . 



The parallel structure of the spatial delta functions and 
the number of viscosities are both imposed by pure sym- 
metry and cannot be disputed. But we may ask which 
of the surface viscosities will have an effect on measur- 
able quantities, such as on the flow profile in nanopores 
or on the motion of a Brownian particle. If we regard 
only forces exerted on the rigid walls — and their correla- 
tion functions — , there are only two viscosities which may 
play a role, namely those possessing a normal direction in 
both contributing stresses. They are easily identified in 
the scheme in Fig. [T] to be rji and 775 , the latter of which 
takes over the role of the bulk shear viscosity, and the 
first one the role of the bulk volume viscosity. It is in- 
teresting to see that the value of 775 is entirely indepen- 
dent on the shear viscosity in bulk. It depends more on 
the hydrophobic/hydrophilic character of the walls. The 
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Figure 10: (Colour on-line) Temporal correlations of the random stress tensor near a wall. Legend and inset as in Fig.JTp. 



surface-equivalent of the volume viscosity, 771, is in both 
cases nearly two orders of magnitude larger than in bulk. 

The finding that the surface viscosities do not have 
values comparable to the bulk viscosities brings us now 
to more profound questions on the dynamical theory in- 
volved: For the Brownian particle, we argue that only the 
first surrounding liquid layers do actually interact with the 
particle. There, the relevant dissipation coefficients are not 
the bulk viscosities but the two viscosities 771 and 775. It 
appears thus as a reasonable to ask why we do not see 775 
instead of 77 appearing in the Stokes-Einstein formula for 
the diffusion coefficient? We have to leave this question 
open at the moment. A possible answer might go in the 
same direction as the discussion of time correlations and 
spatial ordering effects of Ref . [27| . 

A possible answer is certainly linked to the fact that 
the surface viscosities here are local quantities. They have 
been determined using coarse-graining over very small boxes. 
The usual definition of a "viscosity" , however, implies lim- 
its to large boxes and long time integrals (or the limit 
w — >■ after the limit k — )• in Fourier space and time [31 
p. 186]). In this limit, the five surface viscosities cannot be 
identified anymore. Another implication of the use of lo- 
cal viscosities, which are not viscosities in the usual sense 
of a large-box limit, is that they do not lead to differ- 
ential equations. They rather appear as integral kernels, 
nonlocal both in space and in time, in integro-differential 
equation. We cannot say much about the possibility to 
find an to interpret such equations at the moment, the 
less as we expect the long-time tails in the autocorrelation 
function, thus nonlinear terms, to intervene. Neverthe- 
less, we would like to stress that we do not expect the 
qualitative nature of the anisotropy described here — that 
is the large difference of surface and bulk viscosity values, 
and that the surface viscosities do not depend on the bulk 
viscosities — to depend on the precise temporal behaviour. 
The long-time tails appearing in the correlation functions 
should only change the quantitative values of the viscosi- 
ties, leaving the presented material untouched. 
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